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Abstract 

We propose a test of Monte Carlo Parton Cascade models based on 
analytic solutions of covariant kinetic theory for longitudinally boost 
and transverse translation invariant boundary conditions. We com- 
pute the evolution of the transverse energy per unit rapidity for typ- 
ical mini-jet initial conditions expected in ultra-relativistic nuclear 
collisions. The kinetic theory solutions under these conditions test 
the models severely because they deviate strongly from free-streaming 
and also from ideal Euler and dissipative Navier-Stokes hydrodynam- 
ical approximations. We show that the newly formulated ZPC model 
passes this test. In addition, we show that the initial mini-jet density 
would need to be approximately four times higher than estimated with 
the HIJING generator in central Au+Au collisions at cm. energies 
200 AGeV in order that parton cascade dynamics can be approximated 
by Navier-Stokes hydrodynamics. 
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1 Introduction 



Monte Carlo parton cascade models |T|, 0, ^, ^, ^, ^ are currently 
being developed to calculate multiparticle production and evolution in ultra- 
relativistic nuclear collisions. These are being written in a new Open Stan- 
dard Codes and Routines (OSCAR) format ||T0[| to enable rigorous testing of 
essential components of the algorithms and insure reproducibility of the nu- 
merical results. Given the present large uncertainties in the formulation of 
QCD kinetic theory and in the validity of algorithms employed to simulate 
off-shell propagation, screening effects, color coherence, parton production 
and absorption, and hadronization, it is important to develop and subject 
codes to standardized tests on problems where analytic solutions are known. 
One reason for concern is that cascade models necessarily violate Lorentz 
invariance through the scattering prescriptions involving action at a distance 



n|, 0, At RHIC and LHC energies, the phase space of par- 

tons span over ten units of rapidity and extreme Lorentz contraction effects 
arise which amplify non-causal numerical artifacts. Numerical prescriptions 



have been proposed to minimize such Lorentz violations effects[]T3[. However, 



without further tests it is impossible to ascertain whether the final numerical 
results are in fact physically correct. Different Monte- Carlo event generators 
can give different predictions even for basic observables such as the transverse 
energy per unit rapidity, dE±/dy. 

In this paper we propose a simple test of the multiple collision algorithm 
of such models through the computation of the evolution of dE±/dy in the 
special case of longitudinal boost and transverse translation invariant thermal 
initial conditions. This observable is of interest because it is the simplest 
one sensitive to collective hydrodynamic phenomena. It is known p!6|] that 
for such boundary conditions, Euler hydrodynamics predicts that the final 
observed transverse energy is reduced by a factor of about two relative to 



the initial value due to pdV work associated with longitudinal expansion |T7 
This very basic feature of ultra-relativistic nuclear reactions arises even if the 
initial conditions are highly inhomogeneous and turbulent as shown recently 
in [0. On the other hand, pQCD estimates suggest that the viscosity of a 
quark gluon plasma is very high, and Navier-Stokes hydrodynamic equations 
predict about a factor of two less work is done by the plasma. The observable 
consequence of such dissipative dynamics is the reduction of the transverse 
energy loss predicted by Euler hydrodynamics [|^, |20|. In fact, the parton 
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mean free paths may be so long that even the Navier-Stokes approximation. 
It appears essential therefore to employ a full microscopic transport theory 
to evaluate precisely the magnitude of that transverse energy loss. The newly 
developed parton cascade models provide an approximate numerical solution 
to the underlying covariant transport theory and are thus ideally suited to 
address this problem. 



In P2i |23[| analytic solutions for the evolution of the energy density 



were found using relativistic Boltzmann kinetic theory p4[. In this paper we 



extend those solutions to initial conditions more relevant to the RHIC energy 
domain and propose these to test parton cascade models. 

The discussion is organized as follows: In section 2, we review essential 
elements of longitudinally boost invariant kinetic theory. The free streaming 
solution is discussed for both the ideal inside-outside and local thermal ini- 
tial conditions. The integral equation for the phase space distribution in the 
relaxation time approximation is presented. In section 3, the integral equa- 
tion for the energy density is derived together with its Navier-Stokes limit 
for both constant and scaling relaxation times. In section 4, the integral 
equation for the transverse energy evolution is derived and solved numeri- 
cally. The results show that for the expected initial conditions based on the 
HIJING event generator 0], kinetic theory leads to much less collective trans- 
verse energy loss than Euler hydrodynamics but more than Navier-Stokes. 



We show that at least one of the newly developed models, ZPC|jT5|, is able 
to pass this test. Finally, we show that the ZPC also reproduces the kinetic 
theory solution under more extreme conditions where the kinetic theory so- 
lution approaches the Navier-Stokes solution. As far as we know, this is 
the first demonstration of a parton cascade solution approaching an analytic 
Navier-Stokes result. Physically, this result is of interest because it shows 
that in order to approach the Navier-Stokes regime initial parton densities 
must be at least a factor of four greater than predicted by the HIJING event 
generator for central Au + Au collisions at RHIC. 

2 Boost Invariant Relativistic Transport 

The covariant kinetic equation for the invariant on-shell phase space density, 
f{x,p) is 

{pdx)f{x,p) = {pu){C{x,p) + S{x,p)) (1) 
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where C is the Boltzmann colhsion term, S is the source term, and is 
the collective flow velocity field. This equation greatly simplifies in the case 
that longitudinal boost invariance and transverse translation invariance is 
assumed pT|, |22| , p3| . In that case, the on-shell phase space distribution func- 
tion depends only on the reduced phase space variables (t, ^,px), where 
t"^ = — and ^ = ri — y'm. terms of the kinetic rapidity y = tanh~^(p^/p'') 
and pseudo-rapidity r] = tanh~^(z/t). The collective fiow velocity field in 
this case is u'^ = x^/t = {chi], 0±, shrj) with (pu) = m_Lch^ and the transport 
equation reduces to 

ipd.)fix,p) = mx(^ch(0|:-^^)/(r,e,Px) 

= M^/(^,e»,Px) • (2) 

The characteristic function ^(r) = ^{t;tq,^) satisfies 

d^/dr = -r"4anh^ , (3) 

and passes through ^ = ^{tq; tq,^) at r = tq. This trajectory satisfies 

tshe(t;r,0 = rshe 
che(t;r,0 = (l + iT/t)hh'd'^' . (4) 



We will drop the r or ^ arguments on ^{t) when confusion cannot arise. 
In the source of partons was modeled by a Schwinger pair production. 



That mechanism models soft beam jet fragmentation for lab energies < 200 
AGeV. At ultra-relativistic collider energies y/s > 100 AGeV, perturbative 
QCD minijet production is currently thought to dominate the parton pro- 
duction mechanism 0, |20[ . Given an initial distribution of massless partons, 
e.g. from HIJINGi, 

9{y, P±) = dNg/dyd^pi_ , (5) 
the first problem is to construct the source phase space distribution, S{x,p). 
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2.1 Free Streaming Case 

To construct the source, consider an ensemble of on-shell partons with pro- 
duction phase space coordinates, {x^,p^]p1 = rn?}. The general form of the 
covariant free streaming phase space distribution is given by [^] 

where P{s,p) is introduced as the Lorentz scalar probability that a parton 
with four momentum p has been formed by proper time t = s m. This 
of course depends on details of the formation dynamics and external four 
vectors and parameters specifying the reaction. 
The on-shell constraint leads to 

6\p-pa) = 2eipo)6ip'-m')5iy-ya)5\p^-p^a) , (7) 

and thus the reduced (6+1) dimensional phase space density, /, is defined 
through 

J^ix,p) = 2e{p^)5{p^ - m^)f{x,p) . (8) 
The on-shell phase space source function in eq.(P is then given by 

S{x,p) = -^{pd^)f{x,p) . (9) 
[pu) 

For boost invariant boundary conditions, 

S{r,^,P±) = ^/(r,e(r,r',0),P±)|r'=. • (10) 
The local energy-momentum tensor and parton current are obtained as 

m 



T^'^ix) = J d^pp^p^Tix^p) = J ^p^p^fix^p) 
r{x) = Jd'pp^J^{x,p) = J^p^f{x,p) . (11) 



The proper energy density and number density are in turn given by e{x) = 
u^UyT^^^{x), p{x) = u^J^{x) We use here the Landau definition of the 
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normalized fluid four velocity u^^{x) = T^^u^/ {uTu)^ such that the energy 
flux, T*^*, vanishes in the comoving frame where = (1,0,0,0). 

The central observable of interest is the parton invariant momentum dis- 
tribution as given by the Cooper-Prye formula 

E-—^ / dE,p'^f{x,p) (12) 



Here T,'^ {(1,(2,(3) specifies the 3D hypersurface on which the momentum 
distribution is measured. The volume element is given by 

d^^ = e^^^{dEyd(,){dE^/d(2){dEyd(s)d(id(2d(s . 

For a fixed time measurement, we take (Ci, C2, Cs) — (x±,-z), E'' = x'* = 
{t, z,x±) so that dE/j, — {d^x, 0). The inclusive parton distribution at time t 
is then given by 

and the transverse energy per unit rapidity at time t is 

dE±{t) f ,2 dN 



dy J dyd'^p 

On the other hand, for the physically more relevant longitudinal boost 
invariant boundary conditions, the final parton distribution requires an inte- 
gration over a fixed proper time freeze-out surface. For r — {t'^ — z^Y^'^ — t/, 
S is conveniently parameterized by = (rjclir/. rjsh?], x_|_) with dTJ^ — 
[Tfchr], —Tfshr], 0)dr]dx±. Therefore, dTi^p^j^ — {Tf'm±ch$,)dr]d:x.±. In this case. 



dN{Tf) 



— j Tfdr]d^x± m_Lch^ f{x,p) , (15) 



dyd^p 

and the basic calorimetric observable is 

e±(r/) ^ = / TfdridW J d'p^mlch^ f(x,p) (16) 

Our main focus here is on how e_L(r) evolves to its final form given an initial 
boost invariant value at r = Tq. 
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2.2 Ideal Inside-Outside Correlation 



In nuclear collisions at ultra-relativistic energies, the initial parton coordi- 
nates are approximately localized on a transverse sheet, with = (0, x_L(j, 0) 
with x_La distributed over a transverse area according to a density p(x_l) ^ 
9{R — r)/nR'^. Transforming to the (r, ^ = r] — y, p±) coordinates 

6^{x-Xa-ps) = -^^5(s - r/m_L)5^(x_L - X_La - — 7") • (17) 

m_|_r m_|_ 

Neglecting transverse expansion, the approximate ideal inside-outside corre- 
lated distribution is given by 

/o(r,e,P±) ^ -^^P{rm^)9{p±) (18) 

The influence of the formation physics can be studied by comparing two 
analytic models of P = P{m^T): 

P(m , r) - / - ^^^^^ ^ fl9) 

^ \ (mxr)V(l + (mxr)2) Model II ^ ^ 

Model I is simplest for analytic tests, while Model II is more realistic based 
on a study ||25| of induced radiation in multiple collisions. This is obviously 
one of the most uncertain aspects of parton cascade models since formation 
physics is beyond the scope of ballistic kinetic theory. At the least, kinetic 
theory could be supplemented by classical fleld evolution such as proposed 
in|]26[|. (See for example for solutions to coupled fleld and kinetic equa- 



tions in the chro mo-hydrodynamics limit.) 

The phase space source function given (|^) is 

d Siar))Pim^T) ^(OP(mxr) 
So{x,p) = - 9{P±) = ^2 ^(P^) ' (20) 

where we used 

smr,o) = -m ■ (21) 

r 

For this ideal ^ = correlated case, the energy-momentum tensor has the 
form Tq" {x)u^u^ eQij-) + 5^ Pt{t) where = diag(0, 0, 1, 1). The proper 
energy density evolves as 

.„(r) = (uTu) = . (22) 
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where the brackets are used to indicate the transverse momentum integration 
(■■■)= / (fp± g{p±) ■ ■ ■ Note that for the ideal correlated case, the proper 
longitudinal pressure vanishes, while the transverse pressure, Pt = e/2, is 
50% larger than in thermal equilibrium for massless partons. 
The hydrodynamic equations follow from {{dT)u) = a where 



o- r 



j di(fp^mldieS{T,^,p^) (23) 



is the proper energy density source per unit time. The proper energy density 
evolves as 

^ + I = ,M = M%I)> (24) 

The well known solution is given by e(r) = eoij) from eq. (|22|) . For the 
Lorentzian (model II) formation probability in eq.(|l9D, e increases initially 
linearly with proper time until tq ~ 1/ {m±) and then eventually decrease as 
1/r due to one dimensional longitudinal expansion. 
The inclusive distribution (p!3D evolves in this case as 



dyd?pi_ 



ch(y) j dz^-^^^^—y^P{m^r)g{p^) 

showing that partons are formed according to a time dilated formation prob- 
ability. The transverse energy (^1|) develops as 

= (m^P(m^t/ch(y))) (26) 

On the other hand, along a freeze-out proper time surface t = Tf, the inclu- 
sive distribution (^) is given by 

Thus g is the asymptotic invariant parton distribution in the free streaming 
case. The final transverse energy per unit rapidity (|1^) is 

dE (t ^ f f 
e±(T-/) = — ^^-^ = y r/dr/d^xxy rf^Pxm^ch^ /(x,p) 

= {m^Pim^Tf)) = iynR') eoiTf) (28) 
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On account of (p^, dE±/dy tends toward a constant after the formation 
era. The last hne is just the famihar Bjorken formula relating the observable 
dE±/dy to the energy density at the freeze-out proper time Tf. 



2.2.1 Free Streaming with Thermal Initial Conditions 

In order to study the influence of more realistic finite width of the rj — y initial 
state correlations, it is instructive to start with local thermalized source 

Stdr,^,Pi.) = Sir - ro)F (^^^ (29) 

where, eg., F{x) = ce~^. For such a source, the solution of the free streaming 
transport equation is 

/i(r,e>;ro,0,P±) = ^(r-ro)F(^^^ , (30) 

and therefore inverting with (^), 

/i(r,e,P±)=^^(r-ro)F(^che(ro;r,o) , (31) 

where 

che(ro;r,0 = (l + (r/ro)We)'^' . (32) 
The free streaming energy density in this case is 

ei(r) = /rfet^'p±plch2e/i(r,e,Px)=^(r-ro)e(ro)^/i(ro/r) ,(33) 

where the function h modulating the Bjorken tq/t factor is the same as the 



one defined in 22, 23 



°° dC, ch^^ X f sin ^ y/1 — x'^ 



xh{x) = w{x) = X / TTT^ 2 ~ o + 



oo 2 (x2 + sh^O^ 2 V J 
l-^{l-x) + ^-{l-xf + --- . (34) 



For large times, x — >■ and 



TT 

xh{x) —X . (35) 
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The factor 7r/4 is characteristic of the difference between the isotropic and 
inside-outside correlated phase space. The energy density at late times fol- 
lows Bjorken formula reduced by this factor, e(r) ~ f e(To)(To/''")- The initial 
local thermal energy density is given by 

e(ro) = J d^d^p^ pldi^^F{p^cH/To) = KsbT^ , (36) 

where the Stefan-Boltzmann constant is 

POD 

KsB = 47r / dxx^F{x) . (37) 







The transverse energy per unit rapidity is closely related to the energy 
density as seen in eq. ([T6|) . 



ei(r) = {tt^R^) j d'qd^p^plcHHr.^.VA.) . (38) 

This differs from ei(r) in eq.(^) by only one power less of ch^ in the inte- 
grand. Noting that 



ooch^e(t;r,0 4r 

the transverse energy per unit rapidity remains a constant during free stream- 
ing 

ei(r) = ^(ro7ri?2)e(ro)=;im(r7ri?2)ei(r) , (40) 

unlike the energy density. The Bjorken relation between e^ and e is recovered 
in this case only at large times, after the initial local isotropic conditions 
evolve toward the ideal ^ = case. 

2.3 Collisions in the Relaxation Time Approximation 

The simplest way to extend the free streaming analysis to include the effects 
of collisions in eq.(|ID is via a relaxation time approximation pT|, p2| : 



C{x,p) = -{pu){f-h,)/T,{x) (41) 
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where the equihbrium distribution, feq{x,p) is constrained by energy mo- 
mentum conservation {dTu = 0) to obey 



[^P'ipu)if-fe,)/n = . (42) 



po " 

While the relaxation time approximation is strictly valid only for small de- 
viations from local equilibrium, it provides one important limit where the 
general nonlinear transport equations encoded in parton cascade models can 
be subjected to an analytic test. 

With longitudinal boost invariant boundary conditions the equilibrium 
distribution depends only on {pu) and the local temperature T(r), 

U.,,^F[l^)-f„{;^) . ,43, 

Eq. (|4^ ) then constrains the temperature to evolve as 

e(r) = ee,(r) = KsbT\t) (44) 

as long as Tc is assumed to be independent of p. 

The transport equation along the characteristic ^{t; r, ^) that passes through 
^ at proper time r is 

^/(t, e(t; r, 0, P±) = -7c(r)(/ - fe,) (t, r, 0, P±) + S{t, ^(t; r, 0, P±) ■ 

(45) 

This is solved[22] using e ^dr{e^f) to shift 7c/ to the left hand side of the 
equation taking 

/3,(r)= rdtj^it) , (46) 



in terms of the collision rate 

1c{t) = T-\r) = ^{t) (47) 

The collision rate is assumed to vanishes at r = for our boundary condi- 
tions. 

The integral equation for the phase space density is therefore 



(48) 
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Defining the scattering kernel K{t, r) via 

ir(t,r)=/3,(t)e^^W-'^^(^)=7c(t)exp(^-^'-^j , (49) 

the integral equation for / can be written as follows: 

/(r,e,P±)= dtK{t,T) (f,,{t,at;T),p^) +T,{t)S{t,at;r),p^)) . (50) 

Note that the kernel tends toward 6(t — r) in the rapid thermalization limit, 
Tc(t) — > 0, and that in that case / — >• f^g in the source free region. 

For the special case of boost invariant thermal initial conditions, eg. (|29|) , 
/ evolves for r > Tq according to 

/(r,e,P±) = ^(r-r„)e^^(-°)-'3=WF(p^che(ro;r)/To) (51) 
+ r dtK{t,T)f,g{p^chat;T)/T{t)) . 

J TO 

This also shows how the system forgets its initial condition as it tries to 
evolve toward local thermal equilibrium. Whether it get close to feq or how 
far it lags behind depends of course on the functional form of the collision 
rate 7c (i") • 



3 Evolution of the Energy Density 

The integral equation for the proper energy density is obtained by integrating 
eq. (|50|) using 

e(r) = I d^dWpW^ f{r,^,P±) • (52) 

This leads to 



e r 



) = £ dtK{t,T) J d^d'p^plch^^ (^feq{t,^{t-T),p^)+T,{t)S{t,^{t-T),p^ 

(53) 

The transverse momentum and ^ integration over f^g can be done as in 
section2.3 yielding a factor K sBT^{t)w{t / t) . The main simplification arises 
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as a result of the constraint (|^) that allows us to replace KsBT^if) by e(t). 
Therefore, 

e(r) = r dtK{t,T){e{t)w{t/T) + T,{t)cy{t)) , (54) 
Jo 

where the source term depends on the form of the formation ^ correlations 
via 

ait) = J d^d'p^plch'^ Sit, 0, P±) ■ (55) 

For ideal ^ = correlations, a, reduces to eg . (^Sf) . For a thermal correlated 
source as in, eg. (p9D , a reduces to ait) modulated by the same wit/r) factor, 
(HD, as the e(t) term in (H): 

e(r) = r dtKit,T)wit/T)ieit) + T,it)ait)) 
Jo 

= 9ir - ro)e(ro) rc(ro)i^(ro, r) w;(ro/r) 

+ r dtKit,T)wit/r)eit) . (56) 
Jo 

In the free streaming limit, Tc oo, K ^ but TcK ^1. In that limit 
e reduces to ei in eg.(p3D. In the opposite, rapid thermalization limit (^) 
reduces to hydrodynamics as we show in the next section. 



3.1 Navier-Stokes Limit: Constant Tc 

In the limit of rapid thermalization, Tc — > we can expand eg. (|56|) in a power 
series in Tc. Consider first the case where Tc independent of time|2l|, |2^. In 
this case 

K(t,r)=e(*--)/-/r, . (57) 

The kernel is highly peaked at t = r. Thus for slowly varying functions. Fit) 
we can systematically expand 

r dt Kit,T)Fit) = r''"'rfxe-^(F(r)-r,a;F(r) + ^r,VF(r) + --- 
JO Jo 2 

= F(r) - r,F(r) + r,2F(r) + • • • (58) 
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The formal expansion of (|56| ) to second order in Tc then leads to 

e(r) = e(r)-r,(e + |--a(r)) (59) 
9,.. 8 e 4 e . 4cr, 

where we w{l) = l,w{l) = 4/3, w{l) = 4/5 has been used. 
The full integral equation eq. (|56D reduces in this case to 



To lowest order in Tc, (|60D reduces to the Euler hydrodynamic equations for 
Bjorken boundary conditions 

e ^ -^- + a + 0(re) . (61) 
o r 

Thus to lowest order 

The first order correction term can therefore be simplified using 



.. oe 4e 44e . Aa ^, ^ 

^ + Q - + 7- = oT7- + ^ + o - + Or, . 63 
3 r 5 3 15 3 r 



Note that the a terms in (|60|) cancel, and eq. (|60|) reduces to the Navier- 
Stokes, viscous hydrodynamic equation 



^ + \^ = \^2^ + 0{rl) , (64) 



with the expected p^, |T9[ shear viscosity coefficient given by 

4 

V = Y^r.e . (65) 

For constant it is clear that viscous corrections can be neglected at large 
times. 



14 



3.2 Navier-Stokes in the Scaling Limit: rc = r /a 

For the Bjorken boundary conditions, however, the proper density p = (uJ) 
decreases as p(r) = tqp{tq)/t due to longitudinal expansion. Consequently, 
the collision rate is expected to scale[^] as 

7c = = {(^tp{r)) =a/T (66) 

The coupling parameter that controls the rate of thermalization is in this 
case approximately constant 

a = (o-tp(ro)ro) = ro/rc(ro) (67) 

where at is the transport cross section, and tq is any time after the source 
is negligible. For a finite formation probability, such as in eq.(|T9|), a must 
initially grow from zero. In this section we ignore this initial time dependence 
of a to simplify the analytic treatment. In this scaling regime, the kernel is 

Kit,T) = ^('-y (68) 



t \T 

The integral equation (|56| ) for the energy density reduces to 

e(r) = alj (-) h{t/r){e{t) + ta{t)/a) (69) 

For a ^ 1 the kernel is strongly peaked near t = t, and we can expand, 
analogous to eq. ( ^81) , in powers of 1/a 



r dt K{t,T)F{t) = a dxx''-^F{r -t{1 -x)) 
Jo Jo 



F(t) — H — \ 

^ ^ a + 1 (a + l)(a + 2) 



(70) 



The source term expanded to second order in this case is 



'— + 0{l/a') , (71) 



15 



while the energy density term expands to 



At 0(1/q;), we again recover the Euler hydrodynamic equation (pT|). To next 
order we recover the Navier-Stokes equation (0) with a time dependent 
viscosity 

This has the same form as in (^) except that Tc is here time dependent. 

The most interesting point associated with the scahng Navier-Stokes is 
that unhke in the constant Tc case (|6^) , the viscosity term does not become 
neghgible compared to the pressure term ( |e/r) at late times. Dissipation 
in this case decreases the effective speed of sound from = 1/3, appropriate 
for an ideal ultra-relativistic gas, to 

As discussed in P3[, this implies that the energy density in the source free 



region always decreases more slowly than Euler hydrodynamics predicts: 

.M=.(r„)(^)'"^ . (75) 

It is important to emphasize however, that the Navier-Stokes approxima- 
tion can only apply if the viscous term is small compared to the pressure 
term. This requires 

a dN 

This condition is certainly violated at early times when dN/dy is small. In 
addition, as we will see, this condition is not satisfied in a dense parton gas 
if the screened cross sections are as small as pQCD estimates would indicate. 
The detailed study of deviations of kinetic theory and parton cascade solu- 
tions from Navier-Stokes solutions is useful as a gauge of whether particular 
observables such as transverse energy can be correctly interpreted in terms 
of thermodynamic concepts. 
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4 Evolution of e±{T) = dE±/dy 

The transverse energy per unit rapidity is related to the phase space density 
via eq.(^). For the case of local equilibrium, 



J rdi(f^^ J d^p^plchUeqiPLch^/nr)) = ^(rnR')e,^{r) . (77) 



The integral equation for e_L is analogous to (^) involving one power less of 
ch^ in the integrand 

e^(r) = dtK{t,T) J drid^p±{TTiR^)plch^ (78) 

Because of eq.(^) the integral over gives / K{t^T^6,^{t) without the 
weight function w of eq. (|3^) . 

For local thermally correlated sources, in particular, the integral equation 
reduces to 



ex r 



(r - ro)ex(ro) n{To)K{T,,T) + f dtK{t,T)e'l{t) , (79) 



Unfortunately, we cannot replace e^(t) by e^{t), as in the case of the energy 
density evolution equation, because energy conservation only constrains e = 
eeq- Therefore, we must first solve the energy density integral equation (^) 
to compute the temperature, T(r), and then use that solution to compute 
e^(r) via (^). Numerically, it is simplest solve (^,^) together in the same 
loop iterating time steps and keeping track of prior time steps in an updated 
array of e(r„) 

e(r„) = T,{r^)K{r^,T)w{-)t{r^) + Y./\rK{n,Tn)w{-)t{n) (80) 



e±(r„) = -rf |^r,(ro)ir(ro,r)roe(ro) + ^ArK(ri,r„)rie(ri)j .(81) 



This system converges very rapidly even for arbitrary time varying relaxation 
rates, 7c (t), not hmited to the scaling form, (0). Because of its numerical 
simplicity we still refer to the solutions of the above system as analytic for 
tests of parton cascade codes. 
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5 Illustrative Numerical Tests 



In this section we solve the above integral equations for ex(r) and compare 
them to hydrodynamics and parton cascade results for an initial longitu- 
dinally boost invariant but thermally correlated parton gas with e±{To = 
0.2 fm) = 484 GeV and very high parton density p{tq = 0.2 fm) = 20/fm^ 
initial conditions for RHIC energies as estimated from HIJINCp] simulations 
for central Au + Au collisions aX ^/s = 200 AGeV. 

First we show in Fig. 1 the expected hydrodynamic evolution in the case 
of Euler non-dissipative and Navier-Stokes approximations. Note that the 
Euler solution decreases by over a factor of two as expected from previous 
studies. The transverse energy loss is however reduced as the transport cross 
section is reduced. For at < 2mh (corresponding to a < 0.8) the Navier- 
Stokes approximation fails badly even for this very high initial parton den- 
sity. In the figure we have in fact replaced the unphysical negative effective 
pressure in the Navier-Stokes calculation by zero. 

In Fig. 2 we compare the kinetic theory solutions of (|8l|) to Navier-Stokes. 
For cr = 32 mb (a = 12.8) the transport and Navier Stokes solutions are 
almost identical. However, by cr = 3 mb, the Navier-Stokes approximation 
significantly under-predicts the work done during the expansion relative to 
the kinetic theory solution. The kinetic theory result of course also differs 
significantly from the ideal Euler hydrodynamic result. 

In Figs. 3 and 4 the kinetic theory solutions are compared to one of the 
newly developed parton cascade models, ZPC 0. In that model, only elastic 
scattering is currently implemented with a differential cross section of the 
form 

da/dt = + 1)1 [t - /i2)2 (82) 

2 s 

The total cross section is thus independent of energy, a = 97ra;^/2/i^. How- 
ever, the more relevant transport cross section needed for the input to kinetic 
theory (0) is 

= 9^(1 + h) L —{y- )(s + ^^^-y) (83) 

z, s J y 

With this s dependence, the average transport cross section changes with 
time. Numerical results of ZPC for the time evolution of this quantity can 
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be parameterized by 



^ ^ [ 1.80 - 0.255/ mb fig.4 ^ ^ 

For the physically reasonable values of as = OAT and /i = 3 fm~^, the 
transport cross section is so small that even though the parton density is 
very high, the scaling Navier-Stokes parameter a = 0.7 is small and thus 
neither Euler nor Navier-Stokes approximations apply to this problem. The 
essential point is, however, that the ZPC code successfully reproduces the 
analytic results. 

We note that in these examples the ZPC code was used in the transverse 
periodic boundary condition mode to eliminate effects of transverse expansion 
and therefore be directly comparable to the simple analytic results that were 
obtained assuming transverse translation invariance. Numerical results with 
finite transverse radius nuclei, show a reduced transverse energy loss that 
can also be qualitatively understood as due to a more rapid decoupling or 



freeze-out that occurs in that case [53]. With periodic boundary conditions 



the system continues to interact for a longer time. For the conditions shown 
each event takes about one hour to run on a Sparc 20 from tq = 0.2 to r = 6 
fm. The numerical results (solid dots) correspond to averaging 20 events. 

In Fig 4 the results obtained by increasing the initial parton density by 
two (by decreasing the periodic transverse factor of two) are shown. 

These initial conditions involve one half the transport mean free path relative 
to fig. 3 and correspond to the scaling Navier-Stokes parameter a ~ 1.4. For 
the kinetic theory curves the actual time dependent rates of ZPC using (^) 
were used as input in solving (^Tj). The accurate reproduction by ZPC of 
the kinetic theory results in this case confirms again the numerical accuracy 
of the ZPC model in cases far from both free streaming and hydrodynamic 
limits. 

Finally in Fig. 5, we show results where the initial density was again 
doubled to consider with a 3 that should closer to the Navier- 

Stokes regime. In this case we again find that the ZPC and kinetic solutions 
again coincide remarkably well, and in this case the viscous hydrodynamic 
solution approximates much better the kinetic theory evolution. Physically 
the price paid for this agreement is the necessity to increase a oc (Jtp{To) 
by a factor of four relative to the HIJING estimates. This can be achieved 
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physically in many ways. One way is to increase the initial mini-jet density 
by a factor of two (by decreasing the mini-jet scale po in HIJING from 2 
to 1 GeV) and at the same time reduce the screening scale from /i ^ 0.6 
GeV to approximately 0.3 GeV. This however is strongly contrary to pQCD 
expectation where /x ~ gT should be an increasing function of the den- 
sity. In addition such a small po is inconsistent with the data on pp at FNAL 
energies 0. One could expect that inelastic gg ggg--- processes, not yet 
included into ZPC could increase gluon density somewhat. However, for our 
chosen initial conditions, dE±/dy{TQ ~ 0.2 fm) ^ 0.5 TeV in Au + Au with 
p(ro) ~ 20/fm^ in Fig. 3, the initial energy density is e(ro) ~ 17 GeV/fm^, in- 
stant chemical and thermal equilibrium would imply that the density would 
be pth{To) ~ 2Tq ^ 17/fm^. Thus, the initial conditions are in fact close 
to chemical equilibrium in this case and further gluon multiplication is not 
likely. It thus appears difficult to approach the Navier-Stokes regime at RHIC 
energies. At LHC energies initial densities up to a factor of ten higher par- 
ton densities may arise and thus collective hydrodynamic behavior should be 
more easily be achieved. 

Aside from providing physical insight into the possible collective behav- 
ior parton evolution in nuclear collisions, the above examples serve well to 
illustrate how our proposed analytic test can be applied to newly developed 
cascade models as a check of the numerical implementation of kinetic theory. 

6 Summary 

The evolution of the transverse energy per unit rapidity was computed in 
kinetic theory for initial conditions of partons that may arise at RHIC en- 
ergies. This observable is one of the basic probes for collective phenomena 
and is particularly interesting because Euler hydrodynamic equations pre- 
dict a factor of two loss associated with work down as the plasma expands. 
However, with pQCD estimates of the density of mini-jets and the screened 
pQCD cross sections it appears that local thermal equilibrium may be hard 
to maintain even if the initial conditions are assumed to be in boost invariant 
local thermal equilibrium. Kinetic theory provides a microscopic transport 
theory to evolution far from local equilibrium. The newly developed parton 
cascade codes therefore appear to be essential to compute realistic signatures 
in nuclear collisions. 
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In this paper we proposed a simple analytic test of parton cascade models 
for idealized longitudinally boost and transverse translation invariant initial 
conditions, and showed that this test could be implemented at least in the 
ZPCpl code. In the course of applying this test, initial discrepancies between 
the analytic results based on eqs. (|8T|) and numerical ZPC results were found 
to be useful in debugging the code. A similar test|2^ of the preliminary 
version of the VNI codeQ] uncovered an unphysical scattering prescription 
that led apparently to anti-work. In the final version of that code that pre- 
scription was corrected by modifying the low transverse momenta parts of 
the scattering subroutines so that similar results to the ZPC test were found. 
The fact that analytic results far from the Navier-Stokes domain could be 
reproduced by several parton cascade codes is an important step in demon- 
strating the soundness of the general parton cascade numerical technology 
being developed for applications to nuclear collisions at RHIC and LHC en- 



ergies. Our main proposal is to subject all codes at least in the OSCAR|T0 
repository to such tests. 
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Figure 1: Comparison of the evolution of the transverse energy per unit ra- 
pidity in ideal Euler hydrodynamics versus Navier-Stokes dissipative hydro- 
dynamics. The initial conditions were fixed at r = 0.2 fm to be dE±/dy = 482 
GeV, p = dN/dy / {tA±) = 20//m^. The dependence on the transport cross 
section via the viscosity coefficient is shown. 
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Figure 2: Comparison of kinetic theory evolution to Navier Stokes and Euler 
hydrodynamics for several cross sections given for the same initial conditions 
as in fig.l. Note the slow convergence of Navier-Stokes to the kinetic theory 
result as the transport cross section increases. 
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Figure 3: Comparison of analytic kinetic theory results to numerical ZPC 
code[0 results obtained by averaging 20 events. A periodic transverse grid 
of dimensions 10 fm was used. Initially (at r = 0.1 fm), Tq = 500 MeV, in 
an interval —5 < r] < 5, with ^ = 400. The screening mass was assumed to 
be /i = 3 fm"^ with a strong interaction coupling constant as = 0.47. The 
interaction length for the parton cascade was 0.3 fm, and the initial mean 
free path was ~ 0.3 fm. The comparison starts at r = 0.2 fm. The good 
agreement found with the analytic results confirms the validity of the ZPC 
algorithm in an interesting case where the exact results deviate strongly from 
free-streaming and dissipative hydrodynamics. 
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Figure 4: Comparison of ZPC results with analytic kinetic theory for initial 
conditions with twice as large initial parton density as in Fig. 3. Here 
the scaling Navier-Stokes with parameter a = 1.4, still deviations strongly 
from the exact kinetic theory results. Nevertheless, good agreement with the 
numerical ZPC model is again found. 
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Figure 5: Comparison of ZPC results with analytic kinetic theory and 
scahng Navier-Stokes for initial conditions with the parameter a — 2.8. This 
demonstrates the ability of the ZPC cascade model to approach the Navier- 
Stokes dissipativc hydrodynamic domain under extreme initial conditions 
corresponding to four times the default HIJING parton density. 
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